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I. INTRODUCTION 


A. PROBLEM CONTEXT 

Estimating the direction of arrival (DOA) of a signal is 
a problem of great importance in the operation theater of a 
task force at sea. The direction of arrival (DOA) vector may 
be used to distinguish a friendly approaching aircraft from a 
foe. If the aircraft is outside an expected bearing corridor, 
then it is generally thought to be an enemy aircraft. This 
expectation must be verified by other check points during the 
target identification and classification process. Accurate 
determination of the DOA of a signal is a very important issue 
in the operation of weapon systems and should be completed 
with extreme care. 

The main goal of this work is to investigate the accuracy 
and speed of high resolution DOA techniques. The technique of 
Choice will allow the DOA to be solved accurately in real time 
by an on board computer. 

Passive linear equispaced phased array sensors are used to 
generate signals for the DOA determination. Initially, any 
signal received by the sensors is processed to estimate signal 
direction with consecutive samples tracking potentially moving 


targets. 


B. THEORETICAL BACKGROUND 

Consider the case of a linear equispaced array of n 
sensors receiving Signals emitted from m sources. 
Furthermore, let us assume that the signal may be represented 
by a narrowband signal embedded in additive Gaussian white 


noise. The resulting signal received at the array is 


A (1) 


where s, represents all narrowband components emitted by the 
source(s) and n, is the additive noise. 
The output at the q” array element may be expressed [Ref. 


1] as 


ud 2 -1) :d'sin 9, 
Yea A,exp [j-[y c- S229 Ot ipd ge a 


LOT i cR EHE 


where 0, represents the i" signal arrival angle, w, represents 
the center frequency of the narrowband sources, d represents 
the array element spacing, A, is the amplitude of each incoming 
signal, A is the wavelength related to the traveling wave 
movement, 4$, represents the random phase of each incoming 
signal, assumed uniformly distributed over [0,27]. Finally, 
D, lS a zero mean random variable representing the noise at 
the q” array element. The output vector at time t is a n- 


dimensional column vector 


xc n eas ie (3) 
The mode vector m is 


The noise vector is 


n- (5) 


Equations 3-5 may be used to write the received signal vector 


as 
al: 
E mexp(-2myasim (05) / 2) 
Ajexp lj [w t+0,]] +n,, (6) 
i=1 
exp (-277d(n-1) sin(0,)/A)j 
where the vector x, is defined for continuous time. The 


estimated received signal correlation matrix is given by the 


following equation: 








nest 
Å H E H 
RU deca = oam. (7) 
ai M sede nest 


The vector, x,, is the expression defined in Equation (6), for 
each discrete time interval. The quantity nest is the number 
of estimates used to compute the autocorrelation matrix, R,,. 

The autocorrelation matrix, R,,, may be used to extract the 
signal information by using high resolution techniques such 
as: Multiple Signal Classification (MUSIC) [Ref. 2] or 
Minimum Norm [Ref. 3]. These high-resolution techniques use 
the singular vector decomposition (SVD) or the eigenvalue 
decomposition (EVD) of the received signal autocorrelation 
matrix to estimate the DOA information. A new SVD (or EVD) 
decomposition is needed for each update when tracking moving 
sources. Since the SVD (or EVD) =decomposition is 
computationally expensive, good results are not obtainable for 
real time systems. The goal of this work is to find a good 
approximation of the signal and the noise subspaces and to be 
able to use these approximations to find the desired signal 
information. To allow their use in real-time algorithms, it 
is anticipated that these approximations will be a compromise 
between accuracy and processing speed. 

The method to be used is the Rank Revealing QR 
Factorization (RRQR). Assuming that AII=QR is the classical QR 
factorization of A, where Q is orthonormal, R is upper 
triangular, and II is an orthonormal permutation matrix, this 


method provides a matrix of the form 


11 F1 
R= E 
E R,> 


where the norm two of R,, ||. | 








;, is small when A is ill- 
conditioned. This is accomplished by means of a special 
pivoting scheme following the cited QR factorization. From 
the dimension of the block R, of this matrix, its rank may be 
deduced. This method also provides the means to find 
approximations for the signal and noise subspaces of the 
matrix A. It will be shown that this information is contained 
in the orthogonal matrix Q, generated from the RRQR. Details 
about the method follow in Section A of Chapter 2. 

The RRQR will be used to decompose the noise-free 
autocorrelation matrix. This matrix is obtained by 
Subtracting gl from the signal autocorrelation matrix defined 
mer), where o is the Gaussian white noise variance and I is 
the identity matrix. The need for the noise-free 


autocorrelation matrix is due to the special pivoting scheme 


developed by Chan [Ref. 4]. This scheme works for rank 
deficient or ill conditioned matrices. However, the 
autocorrelation matrix, py does not present this 
characteristic. Therefore, no accurate signal subspace 


information can be obtained from the RRQR decomposition of R,. 
To correct this deficiency of the algorithm, information about 


the Signal to Noise Ratio (SNR) is needed. 


It is possible to update the classical QR decomposition of 
the noise-free autocorrelation matrix. The possibility of 
updating the RRQR will also be investigated, as this technique 
allows tracking of moving targets. The rank-one modification 
of a matrix will be used to update the matrices Q and R 
without accessing the updated noise-free autocorrelation 


matrix. 


ET. THE RANK-REVEALING QR (RRQR) FACTORIZATION ALGORITHM 
This chapter introduces the theoretical background 
necessary to understand the RRQR algorithm presented by Chan 
[Ref. 4]. This decomposition will then be used to estimate 


the signal information. 


A. THE RANK-REVEALING QR ALGORITHM 

The QR factorization is a decomposition of a given m by n 
Matrix A into three other matrices, Q, R and II, so that All=0R. 
The matrix, H e R"", is an orthonormal permutation matrix, Q 
EE 15 orthonormal, i.e., Q"Q=I,, and R e C™ is upper 
triangular. The QR decomposition, using a fixed permutation 
II, provides a unique pair of matrices Q and R. [Ref. 4] 

When A is full rank, R is non-singular. However, if Ais 
rank deficient with rank deficiency r, we may choose II so that 
the rank deficiency of R is exhibited in the form of a small 
lower right block R,=0, as shown in Equation (9). Note that 
since Q and II are orthonormal matrices, their singular values 
are equal to one. Therefore, A and R have the same rank. 


The matrix R may be defined as 


ii Rio 
R= 
É R,, 


where R» is r by r [Ref. 4]. 


(9) 


a 





Let us assume that c, is the i" singular value of A, where 


O,20,2...20 It is possible to show that GPs Ro [Ref. 


nº 





4]. Thus, if the norm two of the matrix Ry, |R»l,, is small, 
A has at least r small singular values. The converse is not 
true, that is, it is not guaranteed that, if A has r small 


singular values, its QR factorization will yield a small 





[Rollo - 
Thus, it is necessary to find an algorithm that is able to 
reveal the matrix rank via a small R, block. The Rank- 
Revealing OR (RROR) factorization algorithm (Ref. 4] provides 
such a decomposition of A. 

In this section, a brief discussion and the theoretical 
background of the method is provided. The main purpose here 
is to identify the characteristics which are important to the 
Direction of Arrival (DOA) problem. A thorough theoretical 
treatment of the RRQR factorization has been written by Chan 
(Ref. 4]. 

This proof starts from a rank-one deficient matrix. It 
will be extended later for a rank-r deficient matrix. Let us 
assume that a given matrix A is nearly rank-one deficient. 
The Theorem 2.1 of [Ref. 4] is reproduced here, adapted to the 
complex case. 

THEOREM 2.1 Suppose that we have a vector x ec C? with 
|x||,=2 such that Jarl ,=e, and let II be a permutation matrix 


such that if II'x=y, then |y,|-lyl.. Then if AII-QR is the 
OR factorization of AI then 


|n TEN 


Proof: First we note that since Iv.l - lvl and ly l2 11x] o1, 
we have 


1 
Va > VE (10) 


Next, we have 


QFAx - QFAIIII7x = (11) 


l 

I 
= 
| 


Inn) 


Therefore, 


eU Axl, = lo “Ax], (12) 





IRyt> [z Y us 


f 


and from (10) and (12) we have the desired relationship. 
Now, assume that v, e C^ with |v,|;1 is the smallest right 


singular vector of A, then we have 


|Av]|., =0,,. (13) 


If we define 


| (IT7v) |. 2|] vl., (14) 


where II is the permutation as defined in Theorem 2.1, we see 
that AHI has a QR factorization with pivot r, and that the 
magnitude of r, is less or equal to |o; "nl. In other words, 
the (n,n)* element of R is small. Therefore, it is possible 
to make r, small with an adequate choice of the permutation 


matrix. 


Since we need only an adequate permutation matrix Il, an 
approximation of the SVD of A for the smallest singular 
vector, v, is adequate. This point shows the advantage of the 


RRQR algorithm over the SVD. An approximation of v may be 


computed much faster than the true singular vector. In 
Chapter 3, methods to find a good approximation for v is 
demonstrated. 


The above results serve as the foundation to compute any 
QR factorization of A. An approximation of the smallest right 
singular vector of A is found. Then we determine II, as in 
(14), and compute the QR factorization of All. Alternatively, 
the eigenvector corresponding to the eigenvalue closer to zero 
(Ref. 5] may be found, rather than the smallest right singular 
vector v. The proof that this algorithm is valid for the 
eigenvector case is trivial, as c is replaced by [Al 
throughout the proof. 

The above one-dimensional algorithm may be extended to the 
case when A is nearly rank-r deficient, with r>1. We want to 


find a permutation II, such that 





R 
All = oR = o|” 12 (15) 
Ur 
is the OR factorization of A. The submatrix Rk. às (Fa 


dimensional and |R,l, is small. We apply the one dimensional 


algorithm presented above to R,, for r=1,2, 225, where NES 


TO 


the leading principal (n-r x n-r) dimensional submatrix of R. 
After isolating a (r x r) dimensional R, block, we may use 
this one-dimensional algorithm to compute a permutation P such 
that R,P=Q;R, is the QR factorization of R,P. Next, we 
isolate a (r+1 x r+1) block. This guarantees the (n-r,n-r)* 


element of R is small, leading to: 








2 T 
ad Qi Ri» (16) 
On RS 
where 
i-u” % (17) 
0 
and 
=Q E 2 (18) 
0 





We notice that Equation (16) is the QR factorization of Af. 
The complete algorithm is summarized below in steps 0 to 9. 
An implementation of this subroutine using the MATLAB™ 
software may be found in Appendix A. 


0. Compute a first QR decomposition of the noise-free 
autocorrelation matrix. 


"For 1=n,n-1,...,n-r+1, do: 
Mer, be the leading i X i block of R. 


2. Compute the singular vector v e C! of R,,, corresponding 
to the minimum singular value o,4,(R,) with |vl=1. 


TRE 


3. Compute a permutation P ¢ C” such that | (P#v),|=|/P#v].. 
(This means find the maximum absolute value element of the 
smallest right singular vector and swap it with the i” 
element of the same vector). 


4. Assign v= 





é CP Co Che column CÊ PNR 


Å T PO 
5. Let W=P"W, where P= : " 


o Compute QÑ, the QR factorization of RyP. 


a 9 
ME 


^ T 
Ri, Or Ry, 
0 Roo 


7. Let IIJIP. 


ger == 





9. Let R= 








For this algorithm to provide the desired R» block of R 
small in norm, we must have a R, block at step two with a 
small singular value to insure that the (i,i)? element of ÅR, 
is small. At step nine, the (n-i)" (1ast) row of QR, must 
be small in norm. If these two assertions are true, then the 
lower (n-i-«1 x n-i-«1) block R4 in step 9 is small in norm and 
the desired QR factorization exists. To prove the first 
assertion, we reproduce the Lemma 3.1 of [Ref. 4]. Chan's 
derivation for the case of a real valued matrix is adapted 


here for the complex matrix case. 


22 


"Lemma 3.1: Let B e C™ be a matrix containing any subset 
OR columns of A. Then 


a CS O So 
The submatrix R, is the subset of QHAII, composed of its 
first i columns. Q" and II are orthonormal and therefore their 
singular values are one. Using this observation and the Lemma 
3.1 of [Ref. 4], we have 
Onin (Ry) < 09¡(Q0"AM) =0,(A) . 
The above inequality guarantees that R,, has a small singular 
Elie if 0.(A) is small. 
To prove that R, is small in norm, we reproduce here the 


Lemma 3.2 and Theorems 3.1 and 3.2 of [Ref. 4]. 


emma) 3.2: The matrix W =[w.«j ..., Ww] € C computed by 
the algorithm RRQR satisfies the following properties: For 
l=n,n-1,...,n-r+1, 


1) Iwill=1, 

2) (w);=0 for j>i, 

3) Iwil=|w. = 1Nå, 

4) |Alwl|,=ö,< 0,(A), where (5=9mim(Ru) )i- 


1 1 


Theorem 3.1: Let the matrix W € C as computed by the 
algorithm RRQR be partitioned as W=(w,",W,"), where WM e 
C" is upper triangular and non-singular. Then the QR 
factorization of AII as computed by the algorithm RRQR 
satisfies: 


Raul: = Guri |W VE. 
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Proof: Denote the columns of W by (w,,,,...,w,). Define 
YasQ"AIIW € C”=", From property (4) of Lemma 3.2, we have 


Aarti 
MP S elle = y 5 7 S / (Cr) DE < pros (23) 
1=n 


where 0, = ("mn (Ru) ); and |. ||, denotes the Frobenius norm. Next 
the matrix Y can be expressed as: 





WER IV 
Y = Ọ TANW = ro DNO (24) 
K55W, 
which leads to 
Ra | 
[YI,>11%,,W, 1,22 . (25) 
|W" lle 
Combining Equations (23), (24) and (25), we get 
R 
IR> ll D 


o 2 -= 
mro natal / 
I, Lyr 


from which the desired result follows. 


Theorem 3.2: The algorithm RRQR computes a permutation 
II and a QR factorization of A, given by All=QR where the 


elements of the lower (r x r) upper triangular block of R 
Satisfy 


Jel 
e oito, pong e 
= 


< 27 46,/n for n-r<isj<n 


(27) 


Proof: Using Lemma 3.2, we have, for n-r«isj s n 
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j 
o ,2|Allw,|,=|Rw,|,2| (Rw,) ill AG 
=i 


Isolating the k=j term in the sum and rearranging terms, we 
get 


ET 


zi (Ww 4| s e5*5, [ric (w)) gl 
k=i 
From Lemma 3.2, | (w,),|=]](w) |]. = 1/73, we have 


j-1 
kar S o J-9. ele 
kei 


Solving this recurrence in the index j, we get the bound given 
in the first inequality in Equation 27. Using the bound ovk 
< oVn in each term of the sum in the desired result, we get 
the second bound.N 

Equation (27) in Theorem 3.2 shows the bounds for the 
elements of the matrix R generated by the RRQR algorithm. The 
second inequality states the bounds for the elements of the 
block Ra. The factor 2? indicates that one element in R, 
increases with its distance from the main diagonal.  There- 
fore, for large values of the rank deficiency r, the bounds 
may be quite large, as it grows exponentially. Numerical 
simulation cases will be shown later to demonstrate that these 
bounds are overconservative. 

This algorithm was originally proposed to be iterated from 
n until n-r+1, where n is the number of sensors in the array 


and r is the noise-free autocorrelation matrix rank deficien- 


cy. Alternatively, Prasad and Chandna [Ref. 5] proposed to 


DS 


iterate the recursion until n-m+l, where m is the number of 
sources. A faster algorithm is the result of this technique 
as generally there are fewer sources than sensors. The 
drawback of this approach, however, is a loss of precision in 
the estimated signal information, as the algorithm might not 
capture all the rank deficiency of the matrix to be decom- 


posed. This problem will be investigated in Chapter 4. 


B. USING THE RRQR ALGORITHM TO FIND THE NOISE AND SIGNAL 
SUBSPACES OF THE NOISE-FREE AUTOCORRELATION MATRIX 
Assume that we have a linear equispaced array composed of 

n sensors and m sources, with n » m. The theoretical noise 

free autocorrelation matrix, (R,), is r (=n-m) rank deficient. 

However, the practical noise-free autocorrelation matrix is 

nearly rank deficient. 

It is possible to find the signal and noise subspaces of 
an incoming signal using the RRQR algorithm shown in Section 
A. After applying the complete RRQR algorithm, the matrix R 
reveals the rank deficiency of R,. Consequently, the norm of 
R» is small compared to the norm of the rest of the matrix R. 
Note that R is upper triangular. Here, R, has dimension (n- 
r(=m) x n-r), Rə has dimension ((n-r) x r), and R, has 
dimension (r x r). The signal subspace is contained in the 
first m columns, and the noise subspace is contained in the r 


last columns of the matrix Q. The matrix, W, of the RRQR 
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algorithm is not used here as an approximation of the noise 
subspace, as originally proposed by Chan [Ref. 4]. The use of 
this matrix to approximate the noise subspace yielded poor 


results, as shown by Fargues [Ref. 3]. 


C. USING THE GIVENS ROTATION TO REFACTOR R,,P 

A new QR factorization of the R,,P term is presented in 
step six Of the RRQR algorithm. In step one, it is reasonable 
to run a complete QR algorithm for the first QR factorization. 
But as R, is already upper triangular and P is only a permuta- 
tion that changes the position of two rows, we can use 
inclusive Givens Rotations to zero out the few non-zero 
elements of R,,P below its main diagonal. These rotations 
yield a more efficient algorithm to deal with the special 
structure of the problem. 


The complex Givens Rotation matrix is defined as 


ra : (31) 
-conj (s) i 


where c is real, s is complex and such that |c|?+|s|7=1. 





The Givens Rotation matrix is defined such that the 


following equality is satisfied [Ref. 6]: 


c- 


Suppose we have a matrix, A, and we want to zero out a 


B 
o 


(32) 








given element a(k,i), below the main diagonal, i.e., k>i. If 


ny 


we multiply this matrix by the matrix J defined below, we will 
make the element a(k,i) equal to zero. If we repeat this 
procedure for all elements of A below the main diagonal and 
different from zero, the matrix A is transformed into an upper 


diagonal matrix. This is the essence of the Givens Rotation. 


J = 0 0 0 0 
= Doni) 0. -Go O k 
| 
Co 0 1 


The subroutine GIVENS1 used to compute the complex Givens 
Rotation matrix J and the subroutine QRGIV used to perform the 
QR factorization using the Givens rotation may be found in 


Appendices B and C. 
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III. METHODOLOGY FOR APPROXIMATING THE SMALLEST RIGHT 
SINGULAR VECTOR 


Chapter 2 presented the RRQR algorithm and indicated that 
a procedure is needed to find reliable approximations of the 
smallest right singular vector of a matrix. The inverse 
iteration method and the Incremental Condition Estimator are 
introduced to address this problem. The inverse iteration 
method is used to find the smallest singular vectors of a 
matrix. It may also be used to find the smallest eigenvector 
in absolute value. The Incremental Condition Estimator is 
used to find a reliable approximation of the smallest singular 


value and the corresponding singular vectors of a matrix. 


A. THE INVERSE ITERATION METHOD TO FIND THE SMALLEST SINGULAR 
VECTOR 
Step two of the RRQR algorithm is a computationally expen- 
sive part of the procedure [Ref. 7]. It is therefore desir- 
able to find an algorithm that finds the minimum singular 
value and its corresponding right singular vector quickly and 
accurately. 
The algorithm implemented is the inverse iteration method. 
Starting with an initial guess, v, (in this case a vector 
composed of ones), the following algorithm is iterated until 


convergence [Ref. 8]. 


L9 


1. Solve AU Eor ua 

2. Let uia mu csl 

3. Solve A vM HN RCRUM 

4. Let Vi Ta Nao. 

Let the converged v, be v,, and compute u,, and o,,, by: 

Aü-v,, uy-ü/lül., e-1/lül. 

The above algorithm computes the right (v,) and left (u,) 
Singular vectors as well the minimum singular value of the 
matrix A. This algorithm yields good results for a small 
number of iterations, as shown by the numerical simulations 
presented below. Three iterations were used here to estimate 
the singular vectors. | 

To evaluate the results achieved with the inverse itera- 
tion algorithm, we examined four different test cases. One to 
three iterations were run on each test case. In each case, 
the right singular vector was compared to the corresponding 
vector generated by the SVD decomposition computed with the 
MATLAB™ software. This procedure was followed for 1000 
different random matrices. 

Comparisons between approximated and computed singular 
vectors were obtained by evaluating the magnitude of the 


projection of the estimate over the true smallest right 


Singular vector. If the two vectors are parallel, the 
projection is maximum and equal to one. If they are perpen- 
dicular, the projection is zero. The projections were 
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distributed among ten bins, ranging from 0 to 1, as shown in 
Tables 1-4. The mean and standard deviation for each one of 
the iteration steps were computed. 

The first test case used a 10x10 dimensional random matrix 
generated using MATLAB™. The second one used a 10x10 random 
matrix, with the lower triangle was imposed as zero. The 
third test case, used an upper triangular matrix R, obtained 
from the QR decomposition without pivoting of a 10x10 random 
matrix. Finally, the fourth test case used an upper triangu- 
lar matrix R obtained from the QR decomposition with pivoting 
of a 10x10 random matrix. These cases were evaluated in terms 
of performance to verify the adequacy of the method within 
different and perhaps more demanding contexts. The test 
results are summarized in Tables 1 through 5. 


Table 1: INNER PRODUCT MAGNITUDE BETWEEN ESTIMATED AND TRUE 
SMALLEST SING. VECTORS, RANDOM MATRIX, 1000 TRIALS. 
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Table 2: INNER PRODUCT MAGNITUDE BETWEEN ESTIMATED AND TRUE 
SMALLEST SING. VECTOR, RANDOM MATRIX WITH LOWER TRIANGLE 
PORTION EQUAL TO ZERO, 1000 TRIALS. 


| Percentage 
of Projec- 

VET On ng 
between 


[.9 and i 98.8 [93.5 [55.5 | 


Standard 0399 0184 071 
Deviation 





Table 3: INNER PRODUCT MAGNITUDE BETWEEN ESTIMATED AND TRUE 
SMALLEST SING. VECTORS, UPPER TRIANGULAR MATRIX RESULTING FROM 
QR FACT. W/O PIVOTING, 1000 TRIALS. 
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Table 4: INNER PRODUCT MAGNITUDE BETWEEN ESTIMATED AND TRUE 
SMALLEST SING. VECTORS, UPPER TRIANGULAR MATRIX RESULTING FROM 
QR FACT. W/ PIVOTING, 1000 TRIALS. 


| Percentage 
of Projec- 
DEPO rang 


"9823 | 9887 


Standard LIO Ar Ue. 350743 
Deviation 





Table 5: COMPARISON BETWEEN THE FOUR CASES EXAMINED IN TABLES 
1 THROUGH 4. 





between .9 and 
O y 
jverage-or | 

Standard Devi- 

jections 









ZS 


Table 5 shows that the best result was obtained from the 
imposed triangular matrix experiment. However, all of the 
experiments show good results. These results justify the use 
of the inverse iteration method as a low cost alternative to 
approximate the smallest right singular vector in step two of 
the RRQR algorithm. The important point is to find a permuta- 
tion so that the element of the least dominant singular vector 
that presents the largest magnitude be positioned at the i" 
row [Ref. 4]. Therefore, it is not necessary to find the 
exact singular vector for the RRQR algorithm. The source code 
implemented to approximate the smallest singular vector by the 


use of inverse iteration method is shown in Appendix J. 


B. THE INVERSE ITERATION METHOD TO FIND THE SMALLEST 
EIGENVECTOR 
Prasad [Ref. 5] states that the RRQR-based algorithm works 


when uSing an EVD (Eigenvector Decomposition) rather than a 


SVD (Singular Vector Decomposition) of the noise-free 
autocorrelation matrix. The validity of this approach is 
investigated here. Theoretically, it is correct to only 


replace o by |A| in the RRQR factorization proof presented in 
Chapter 2. 

The inverse iteration method may be used to find the least 
dominant eigenvector (the one with a magnitude closer to 


zero). Again, we find an initial guess v, , which may be a 
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vector composed solely of ones. tie wEol lowing. alger thm 


should be iterated until convergence of v. 


im SOlve Ay,=v,,, for y,. 


2. Normalize v,-y,/lv,|.. 


Comparing the two algorithms, we note that this inverse 
iteration algorithm uses the same initial two steps as the one 
used for the minimum singular value. To compare the 
performance of these two algorithms, the same test cases were 
run aS in previous section. However, six iterations were run 
to provide a viable comparison. The results are shown in 


Tables 6 through 9. 
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Table 6: INNER PRODUCT MAGNITUDE BETWEEN ESTIMATED AND TRUE 
SMALLEST EIGENVECTORS, RANDOM MATRIX, 1000 TRIALS. 











Ts 
3 SRS; E 
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| 

4.0 6 Pp. 
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Standard .2758 .2491 .2347 -2 124 .2116 .1986 
Deviation 


Table 7 : INNER PRODUCT MAGNITUDE BETWEEN ESTIMATED AND TRUE 
SMALLEST EIGENVECTORS, RANDOM MATRIX WITH LOWER TRIANGLE 
PORTION EQUAL TO ZERO, 1000 TRIALS. 


percentage 
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Table 8: INNER PRODUCT MAGNITUDE BETWEEN ESTIMATED AND TRUE 
SMALLEST EIGENVECTOR, UPPER TRIANGULAR MATRIX RESULTING FROM 
QR FACT. W/O PIVOTING, 1000 TRIALS. 


Percentage of [One TWO Three | | ] 
Projection Iterat jIterat |Iterat |Iterat |[Iterat |Iterat 
ha. ya DNO .8 | 
283 i 1.4 : | 


E 
s 


3 1.4 L2 1.4 





Table 9: INNER PRODUCT MAGNITUDE BETWEEN ESTIMATED AND TRUE 
SMALLEST EIGENVECTORS, UPPER TRIANGULAR MATRIX RESULTING FROM 
QR FACT. W/ PIVOTING, 1000 TRIALS. 


Percentage 
of Projec- 
Prone lying 


[92.7 [94.6 [95.5 — 
S erre gren AS MT DT . 0949 .0803 .0668 


Bane@ard 
Deviation 
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The results found in Tables 1 through 4 ands+6 through 
are summarized in Table 10 below. Again, the case for imposed 
triangular matrices showed the best performance. 


Table 10: COMPARISON BETWEEN THE FOUR CASES CONSIDERED FOR THE 
SVD VERSUS THE EIGENVECTOR INVERSE ITERATION METHOD. 











[________| SVD erv] sw Jero | svo [BIGV] SVD EIS] 


| a b 


Average of | 
projections .98671.8775|.9996|.9774|.9868 .9379.988 719880 


Standard De- 
lation of .0802/.1986/.0071|.10191.0843.176 71 07423 0 
Projections 


The SVD columns show the results of three iterations of 






the inverse iteration method to find the singular vector 
corresponding to the minimum singular value. The columns EIGV 
are the results of six iterations of the inverse iteration 
method to find the least dominant eigenvector. 

Å cursory look at these results show that the eigenvector 
approximation is worse than those for the singular vectors 
approximation. A hypothesis test based on the two samples is 
therefore tested for each case. An assumption is made that 


each sample came from different populations with each group of 
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1000 projections considered to be one sample. The assumptions 


are 


1. The first sample, for the eigenvector case, is a random 
sample from a population with mean yg, and standard 
deviation O0,. 


2. The second sample, for the SV case, is a random sample 
from a population with mean jy, and standard deviation o,. 


3. Both samples are independent of one another. 


4. The samples are large enough to apply the Central Limit 
Theorem. 


The hypothesis test can be described by [Ref. 9] 


Ho: 7H? 
Hi: up 


One hypothesis test will be conducted for each one of the 
four test cases. At a level of significance of 1%, H, will be 


rejected if z s -2.33, where 


avg(x) -avg(y) 


= 
s? så (33) 
E 
n n 


If Hy is false, we may decide that y, is smaller than py. 
In Equation (33), n is the sample size of 1000, avg(x) and 
avg (y) are the averages to be compared and s, and s, are the 
standard deviations computed from the samples. For example, 
the first test case avg(x)=.8775, avg(y)=.9867, s,=.1986 and 
ELEDSOZOtrom Table 10. Applying these values to Equation 


(33), we find z=-16.12. 
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Table 11: RESULTS FOR THE HYPOTHESIS TEST CONDUCTED 
ON THE COMPARISON BETWEEN THE SVD VERSUS THE 
EIGENVECTOR INVERSE ITERATION METHOD. 


Do 1 | 16.12 Reject H => <H 
TD 2 | 6.87 —Reject Hy -» ji « H; 
4 [98 [Cannot Reject | 










As can be seen in Table 11 above, the first average yg, is 
significantly smaller than p, (at 1% of level of significance, 
z<-2.33), for all cases but case four. Hence, we reject Hg 
for the first three cases. These results justify the choice 
of using the singular vector approximated by the inverse 
iteration method for use with the algorithm. Note that there 
is no theoretical reason for not using the smallest eigenvec- 
tor rather than the smallest right singular vector in the RRQR 
algorithm. They span the same subspace. The reasons for 
choosing the singular vector relies solely on the fact that 
the inverse iteration method yields better results. The 
source code implemented to approximate the smallest eigenvec- 


tor by the inverse iteration method is shown in Appendix M. 


C. THE INCREMENTAL CONDITION ESTIMATOR 
A third method tested here to estimate the singular vector 
corresponding to the minimum singular value is the Incremental 


Condition. Estimator TOE) [Ref 10]. Suppose that 
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A=[a,,...,a,] 1s a (mx n) dimensional matrix and let 0, 2...2 
Onn 2 O be its singular values. The minimum singular value of 
A measures how close this matrix is to rank deficiency. The 


@eneaition number becomes 





O 
k, (A) == (34) 
Onin 
Now, suppose we take the QR factorization of A. This 


algorithm is intended to work from a lower triangular matrix 
L. Since R is upper triangular let us define L=R'. There 
will be no loss of generality here because the singular values 
of any matrix and its transpose are the same [Ref. 10]. 

If we have a n dimensional lower triangular matrix L 
generating one row at a time with an approximate singular 
vector x such that Omin (L) =1/ |||], and a new row (v!,y) of L, we 


Ran obtain 


(35) 


"i 





a. 


ve 


such that o,,(L')«1/|yl, without having to access L again [Ref. 


IN. 





Given x such that Lx-d with lall,=1 and aa DEE, let 
us find s=sin($) and c=cos(¢) such that ly ll is maximized, 


where y solves 


L Q 


L'y = 
F 





ʻi = d' (36) 





y E 


I 


The solution for the above equation is 


Sx 
y-7|c-sq«|^5 (37) 
Y 
where a=v'x and |d'l,=|dl,=1. 
Now define 
JL Past (38) 
-~ 1|' 








where B=v"x"x+0"-1. In this case, we have 


s l (39) 


Assuming yz0, the optimal pair (s,c)' is the eigenvector 





1 
lz = eae 


corresponding to the largest eigenvalue \,,, of B [Ref. 10]. 


Also assuming a#0, we may define 


n=-É and p=n+sign(a) n7+1 (40) 


to obtain A,,,=%4+1. Finally the optimal pair (s,c)' is given 


= 
på+1 


After computing the optimal pair (s,c)!, a new approximate 


as 


u 
= 


(41) 











singular vector, y, may be computed as defined in Equation 


(37) and the smallest singular value of L' by 
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EU RE (42) 
MP 





Onin 


Using the above estimator we ran the same cases as for the 
inverse iteration presented in Sections A and B above. The 


results are summarized in Tables 12 through 15. 


Table 12: INNER PRODUCT MAGNITUDE BETWEEN ESTIMATED AND TRUE 
SMALLEST SING. VECTORS, RANDOM MATRIX, 1000 TRIALS. 


Percentage of 
Projection 
lying between 
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Table 13: INNER PRODUCT MAGNITUDE BETWEEN ESTIMATED AND TRUE 
SMALLEST SING. VECTORS, RANDOM MATRIX WITH LOWER TRIANGLE 
PORTION EQUAL TO ZERO, 1000 TRIALS. 


[percentage of| (| of 
per Er uU RECEN | 
plying s | 


[0 and .1 (0.0 | 
Li and .2 (0.0 | 
Ern: 5 p: 
.3 and .4 — 0.0 — 
.5 and .6 0.3 —. 
16 and .7 (0.3 
L7 and .8 (03 | 

8 and .9 El 


Standard De- 
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Table 14: INNER PRODUCT MAGNITUDE BETWEEN ESTIMATED AND TRUE 
SMALLEST SING. VECTORS, UPPER TRIANGULAR MATRIX RESULTING FROM 
QR FACT. W/O PIVOTING, 1000 TRIALS. 


Percentage of] — 
Projection 
lying hos 
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tandard De= 1079 || 
viation 
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Table 15: INNER PRODUCT MAGNITUDE BETWEEN ESTIMATED AND TRUE 
SMALLEST SING. VECTORS, UPPER TRIANGULAR MATRIX RESULTING FROM 
QR FACT. W/ PIVOTING, 1000 TRIALS. 
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Table 16 summarizes the results found for each one of the 
four test cases, comparing the inverse iteration method to 
find the smallest right singular vector and the Incremental 
Condition Estimator. As before, those numbers represent the 
percentage of projections lying between .9 and 1 of the 
corresponding vectors found via the corresponding approxima- 
tion method over the true smallest right singular vector found 


via SVD. 
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Table 16: COMPARISON BETWEEN THE FOUR CASES CONSIDERED FOR THE 
SVD INVERSE ITERATION VERSUS THE INCREMENTAL CONDITION 
ESTIMATOR METHOD. 


from QR 
decomp. w/ 


opa e 


[Average of | == | 
morojections .98671.90151.9996|.9929|.98681.97051.98871.9859 
Standard De- 

lation of .08021.2041|.0071|.04411.08431.10791.0743|.0797 
Projections 





Table 16 shows the magnitude of the projections of the 
estimated smallest right singular vectors obtained using 
inverse iteration (INV. ITER.) method and the ICE, onto the 
smallest singular vector computed via EVD. 

A hypothesis test was again performed on these results. 
Table 17 presents the results obtained for the hypothesis 


test: 


Table 17: RESULTS FOR THE HYPOTHESIS TEST CONDUCTED 
ON THE COMPARISON BETWEEN THE SV INVERSE ITERATION 
AND THE INCREMENTAL CONDITION ESTIMATOR. 


| (INV.ITER.) « m, (ICE) 


= -. 81 Cannot Reject Ho 
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Tables 12 through 17 show that the results using the 
incremental condition estimator are not as good as those for 
the SV inverse iteration. Therefore, the inverse iteration 
used to find the least dominant singular vectors is preferred. 
The source code implemented to approximate the smallest 
Singular value and its corresponding singular vectors via ICE 


is shown in Appendix I. 
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IV. COMPARISON BETWEEN THE PERFORMANCE OF THE RRQR ALGORITHM 
ITERATED FROM DIMENSION n UNTIL n-m+1 AND FROM DIMENSION n 
UNTIL n-r+1 

Prasad [Ref. 5] states that the RRQR-based algorithm may 
be used to estimate the DOA information when the algorithm is 
iterated from n until n-m+1 rather than until n-r-«1 as in 
Chan's work [Ref. 4]. The relative efficiency of the RRQR 
algorithm using both approaches is compared by running 1000 
trials. The scenario is m-2 fixed sources at 30? and 32? with 
n-10 sensors. Therefore, the theoretical  noise-free 
autocorrelation matrix is of size 10 x 10 and has rank two 
(m=2). The near rank deficiency is (r-n-m) 8. Three SNR 
cases are tested (-10, O and 10 dB) for each one of the 
situations, resulting in six simulations. 

In the first situation, n-r-«1 equals three. In the 
second, n-m+1 equals nine. It is intuitively obvious that the 
second one is much faster than the first, as it will be 
iterated only twice, from ten to nine. On the other hand, the 
smaller number of iterations, the less probable that the upper 
triangular matrix R will capture the near rank deficiency of 
the noise-free autocorrelation matrix (R,). 

For each run, the largest principal angle between the 
Signal subspace computed via the eigenvector decomposition and 


the approximated signal subspace computed via the RRQR 


algorithm is used for comparison. The same is done for the 


Sra 


noise subspace. A vector of ones is expected if the subspace 
generated by the RRQR algorithm is parallel to its correspond- 
ing subspace generated by the true eigenvector. In the case 
that the two subspaces are perpendicular, a vector of zeros is 
expected. 

Recall that the cosines of the principal angles between 
two subspaces F and G are defined as the singular values of 
the product DOS [Ref. 7], where the matrices Q, and Qo are the 
orthonormal matrices obtained from F and G. To find the SVD 
decomposition of this product, the inverse cosine of the 
singular values is taken using the largest angle. The largest 
angle represents a measure of the distance between the two 
subspaces. This measure is used for noise and signal subspa- 
ces. 

Tables 18-20 present means and standard deviations 
obtained for the largest principal angle for signal and noise 
subspaces using the RRQR algorithm where the first QR decompo- 
sition, in step 0, is performed with pivoting. This computa- 


tion is done for SNR -10, O and 10 dB. 
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Table 18: COMPARISON BETWEEN THE RRQR ALGORITHM ITERATED FROM 
n THROUGH n-r+1 AND FROM n THROUGH n-m+1 FOR SNR=-10 dB. QR 
DECOMPOSITION IN STEP 0 PERFORMED W/ PIVOTING, 1000 TRIALS. 


Iteration from n=10 Iteration from n=10 
ntil n-r+1=3 ntil n-m+1=9 (Å of 
(POE ETODS= 33 MO > ia 


ELODS= 5 0 AT SE gle in degrees. | 


Table 19: COMPARISON BETWEEN THE RRQR ALGORITHM ITERATED FROM 
n THROUGH n-r+1 AND FROM n THROUGH n-m+1 FOR SNR=Ô0 dB. QR 
DECOMPOSITION IN STEP 0 PERFORMED W/ PIVOTING, 1000 TRIALS. 
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Table 20: COMPARISON BETWEEN THE RRQR ALGORITHM ITERATED FROM 
n THROUGH n-r+1 AND FROM n THROUGH n-m+1 FOR SNR=10 dB. QR 
DECOMPOSITION IN STEP 0 PERFORMED W/ PIVOTING, 1000 TRIALS. 















Iteration from n=10 |Iteration from n=10 
ntil n-r+1=3 ntil n-m+1=9 (# of 
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flops=1,301,732) - gle in degrees. 
gle in degrees. 


Signal [2.24 4.85 
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Noise 2.24 4.85 4.09 
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Tables 18-20 show that the noise and signal subspaces 
yield identical means and standard deviations. This was to be 
expected as the noise and signal subspaces contain the same 
information. 

Table 21 shows the results obtained for the hypothesis 
test performed to compare the two methods of iteration. 
Table 21: RESULTS FOR THE HYPOTHESIS TEST CONDUCTED ON THE 


COMPARISON BETWEEN THE RRQR ALG. ITERATED FROM n THROUGH n-r+1 
AND FROM n THROUGH n-m+1. QR DEC. IN STEP 0 W/ PIVOTING 


in dB 
DO 21.6 — Reject B» -» ji < à. 
[i0 [19.9  Reject B => m<u 










The results show that for SNR 0 and 10 dB, the null 
hypothesis is rejected. Therefore, the larger number of 
iterations yield better results. However, for the case of -10 
dB the results are meaningless due to the small signal to 
noise ratio. They do not lead to detection of the signal 
embedded in the noisy environment. 

Tables 22-24 present means and standard deviations 
obtained for the largest principal angle for signal and noise 
subspaces using the RRQOR algorithm where the first QR decompo- 
sition, in step 0, is performed without pivoting. An hypothe- 
sis test is not necessary in this case due to the clear 
Gifference in results found for the different method of 


iterations. 
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Table 22: COMPARISON BETWEEN THE RRQR ALGORITHM ITERATED FROM 
n THROUGH n-r«1 AND FROM n THROUGH n-m+1 FOR SNR=-10 dB. QR 
DECOMPOSITION IN STEP 0 PERFORMED W/O PIVOTING, 1000 TRIALS. 















[terabtron “fren h=10|Tteératven tromma=so 
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flops= Reus gle in degrees. 
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Table 23: COMPARISON BETWEEN THE RRQR ALGORITHM ITERATED FROM 
n THROUGH n-r+1 AND FROM n THROUGH n-m+1 FOR SNR=0 dB. QR 
DECOMPOSITION IN STEP 0 PERFORMED W/O PIVOTING, 1000 TRIALS. 
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Table 24: COMPARISON BETWEEN THE RRQR ALGORITHM ITERATED FROM 
n THROUGH n-r+1 AND FROM n THROUGH n-m«1 FOR SNR=10 dB. QR 
DECOMPOSITION IN STEP 0 PERFORMED W/O PIVOTING, 1000 TRIALS. 
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Tables 22 to 24 show that better performance is obtained 
when uSing a larger number of iterations. Note that no 
difference in performance is found using a QR decomposition 
with or without pivoting in step 0 of the RROR algorithm in 
the case using more iterations. This result is not true for 
the second case. This is ciear when Tables 18-20 are compared 
to their correspondent Tables 22-24. A hypothesis test is not 
needed to verify this, due to the proximity of the results. 
Therefore, the QR decomposition without pivoting is preferred. 
This method is less computationally intensive when the origi- 
nal algorithm is used, performed with the total number of 
iterations. Note that the option using only two iterations 
does not present the same performance regarding the pivoting. 
Therefore, if one chooses this option, care should be 


exercised when evaluating the pivoting needs. 
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V. COMPARING THE PERFORMANCE OF THE MINIMUM NORM AND MUSIC 
SPECTRAL ESTIMATORS 


This Chapter investigated the computation of estimates of 
the direction of arrival of signals obtained using the RRQR. 
Two high-resolution techniques, the MUSIC (Multiple Signal 
Classification) [Ref. 2] and the Minimum Norm [Ref. 11] are 
evaluated to verify their adequacy when used with the RROR 
algorithm for DOA estimation. 

Rao [Ref. 12] shows that the Mean Square Error (MSE) of 
the Minimum Norm estimator is smaller than the MSE of the 


MUSIC estimator. The MUSIC spectral estimator is based on 


the orthogonality principle. Therefore, the principal 
eigenvectors {v,,v>,...,V,} span the same subspace as the 
signal vectors {e,,e,...,e,}(Ref. 2]. Thus, the signal 


vectors are orthogonal to all vectors in the noise subspace. 
The power density corresponding to the sources DOA information 
is given by: 


il 


S |e^(6) v; (43) 


J=m+1 


Puuszc (0) = 


where v; is the singular vectors of the autocorrelation matrix, 
m is the number of signals, n is the number of sensors or the 


size of R, and e(0) is the mode vector as defined in (4). 
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The angle 0 is varied in fine steps from 0 to 2m”. When it 
corresponds to one of the source DOA angles, e(0) will be 
EUST to e, i=1,2,..,m. Since the sum in the denominator of 
(43) is over m+1 through n, we have the singular vectors 
corresponding to the noise subspace. By the orthogonality 
principle, the product in the denominator of Pyysc will tend 
to zero and Pyygc Will tend to infinity [Ref. 2]. The result 
is a peak at the source DOA angles. 

The Minimum Norm estimator is based on the estimation of 


0,, the source DOA angles, from the eigenstructure of the 


autocorrelation matrix. Suppose we have a vector d so that 
feat d2,...,d,), where n is the number of sensors in the 
array. If this vector has the property that x "d=0, 
i=1,2,...,m, where m is the number of sources present and x 


the data snapshot vector at instant i, then we can find a 
polynomial D(z) as 


DZ ane, (44) 


1=0 


so that the zeros of the polynomial lie at the elements of the 
mode vector corresponding to the source DOA angles. 

The polynomial roots corresponding to the DOA angle 
locations lie on the unit circle for the autocorrelation 
matrix when additive noise is not present. However, the m 


roots of the polynomial D(z), corresponding to the m sources, 
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lie near but not on the actual unit circle when the estimated 
noise free autocorrelation matrix R, is used. 
The Minimum Norm Spectral Estimator presents the following 


advantages: 


1. The estimates of @,, the source DOA angles, are more 
accurate even at relative low SNR values when compared with 
other procedures. 
2 The n-m extraneous zeros of D(z) tend to be uniformly 
distributed within the unit circle and have less tendency 
for false sources. [Ref. 11] 
By less tendency to false sources we mean that the zeros 
of D(z) corresponding to noise are much smaller in magnitude 
than the ones corresponding to sources. The method imposes 


d,, the first element of vector d, to be equal to 1 and 


requires that the quantity Q in Equation (45) be minimum. 
n 
0=3" |d,|? (45) 
k=1 


The effect of this minimization forces the extraneous n-m 
zeros to be uniformly distributed inside the unit circle [Ref. 
EI s 

Next, let E, be the matrix constructed with the signal 


Subspace generated by EVD or SVD. We partition E, as follows: 
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å (46) 


ESF» 
ae 
where g=[e,, e, ...,em]' has the first elements of the signal 


subspace eigenvectors and É, is the matrix E, with the first 
row deleted. It can be shown that in order to satisfy the 


desired minimization and forcing d,=1 [Ref. 11], we have: 


al 


a (47) 
-Eg*/ (1-9 "%g) 








Once the vector d, representing the coefficients of the 
polynomial D(z), is determined via (47), the roots are 
computed. There are m roots corresponding to the sources that 
present a large magnitude compared to the ones corresponding 
to noise. These roots reveal the desired DOA angles in their 
phase angles. 

Results for both, MUSIC and Minimum Norm spectral estima- 
tors are shown in Figures 1-4, for two fixed sources at 30° 
and 32°. It can be seen from these examples that the Minimum 
Norm spectral estimator starts to resolve the two sources for 
a SNR of 5 dB, while MUSIC starts to resolve only for a SNR of 


7 dB. This agrees with the results shown in [Ref. 12]. 
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COMPARISON OF MINORM AND MUSIC SPEC. EST., SNR=!1 
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Figure 1: Comparison between the DOA estimated by MUSIC and 
MINNORM for two standing sources located at 30° and 32° for 
SNR=1 and 2 dB. 
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COMPARISON OF MINORM AND MUSIC SPEC. EST., SNR 3 
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Figure 2: Comparison between the DOA estimated by MUSIC and 


MINNORM for two standing sources located at 30° and 32° for 
SNR=3 and 4 dB. 
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COMPARISON OF MINORM AND MUSIC SPEC. EST., SNR »S 
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Figure 3: Comparison between the DOA estimated by MUSIC and 
MINNORM for two standing sources located at 30° and 32° for 
SNR=5 and 6 dB. 
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Figure 4: Comparison between the DOA estimated by MUSIC and 
MINNORM for two standing sources located at 30° and 32° for 
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MI. USING THE ADAPTIVE RRQR ALGORITHM TO TRACK A MOVING 
SIGNAL 


In a real case, we are interested in being able to detect 
and track a moving source. In this case, the moving source 
information was sampled every t, interval to provide an update 
of the source information. The array of sensors will receive 
the signals from the sources leading to the data vector 
E oo. a] From vector K we compute the 
autocorrelation matrix, R,,, from (7) and form the noise-free 
autocorrelation matrix, R,, by subtracting the noise informa- 
tion, oil. Next we apply the RROR algorithm for each update 
to identify the signal or noise subspaces. Finally, the 
ME mum Norm estimator may be applied which leads to the 
identification of the m source DOA angles from the n-m 
extraneous noise zeros. This algorithm is presented in this 


chapter. 


A. THE ALGORITHM 

In order to track the DOA of a moving signal, a noise-free 
autocorrelation matrix must first be generated. Next, we 
compute a first RROR factorization and find the matrices Q, R 
and II. The noise-free autocorrelation matrix, R,-R,-o?I, may 
then be updated for each one of the next signal positions in 


time, according to Equation (48) below 


SH 


ld 
Ra =R 


+20X “ey XX “ora E. 
where x is a (n x 1) dimensional vector containing the input 
signals at each one of the n array sensors. The update of the 
noise-free autocorrelation matrix is achieved using a moving 
window where the number of snapshots used to compute the 
noise-free autocorrelation matrix is constant. The new vector 
x is incorporated in the autocorrelation matrix information 
for each snapshot update, while the oldest snapshot 
information is discarded. 

A possible approach in updating the information is to find 
the updated noise-free autocorrelation matrix using (48). 
Next we apply a QR decomposition followed by Chan's [Ref. 4] 
Eu ots Chor EET TE a complete RRQR algorithm. 

However, it is possible to update directly the existing QR 
factorization [Ref. 3] for each new time sample. Suppose we 


dd whose 


want to add a rank-one matrix C=x-x" to the matrix R,, 
QR factorization is known as  ORII.. The new matrix 
Ra "=R „+x: x" will have a QR factorization Q,RJI’, where x-x® 
is the desired rank-one modification. 

The rank-one QR factorization update may be computed as 
follows: 

Ro = RQU.x:x - Q(R«QHx-xHDII' - Qn 
Let w-Q":x. Complex Givens rotations can be used to zero 


out all elements of w except for the first component, leading 


ES 
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GR 


ET AE Pe O OT (49) 
If the same Givens rotations are applied to R, it can be shown 
that H=J"...J|/R is upper Hessenberg. Consequently, we have 


heed, ) (R, tw°x )=H+fa 0 ...  0]':x"-H, also upper 


n- 
Hessenberg. {Ref. 7] 

Next, we use complex Givens rotations to compute 
G,E...G,, H,-R,, where R, is upper triangular. Combining all of 


the above, we have the desired new QR factorization 


Ra =R EQ) RI, where 


CO E E (50) 


nols 
The reader should refer to Golub [Ref. 7] for additional 
detail. 

We apply two successive rank-one modifications to update 
the noise-free autocorrelation matrix for the current time 
sample, finding a new set of matrices Q, R anad I. The first 
rank-one modification incorporates the new data vector to the 
autocorrelation matrix. The second accounts for removing the 
old information. Then, we apply Chan's (Ref. 4] pivoting 
scheme to insure a correct estimation of the noise and/or 
Signal subspaces. 

Next, we identify the signal and noise subspaces. The 
first m columns of Q constitute the signal subspace, where m 
is the number of sources present. The last n-m=r columns 


constitute the noise subspace, where n is the number of 
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sensors. Note that it is usually more efficient to use the 
signal subspace, rather than the noise subspace as it is 
smaller. Finally, the Minimum Norm algorithm is applied to 
estimate the source locations. 

A total of m out of the n roots of the polynomial whose 
coefficients form the vector "d" in the MINORM algorithm 
corresponds to the source locations. It is expected that the 
m roots corresponding to the sources lie near the unit circle 
and the remaining ones have smaller magnitudes. Some sort of 
filtering must be applied to separate the m source zeros and 
the remaining n-m extraneous zeros. Filtering may be achieved 
either by sorting the expected range of source angles and/or 
by sorting the magnitude. Thus, the algorithm corresponding 


to the adaptive case becomes [Ref. 3]: 


1) Initialization 


R,=X:XF-wo"I (note that the noise-free autocorrelation 
matrix is unnormalized) where X is defined in Equation 
(7), Chapter 1 and w is the number of snapshots used to 
compute the correlation matrix. 


RRQR factorization of R, (RJI-QR) 

2) Start updating. For k=1 to number of updates do: 
a) R,(k«1)-R,(k) «x(k«1) x" (k«1) -x(k«1-w):xÜP(k«1-w). 
b) Update the above QR factorization applying two succes- 
sive rank-one modifications to R,(k), leading to: 
R,(k+1)=Q,RJI’ (Alternatively we may find a new complete 
RROR decomposition of the updated noise-free 


autocorrelation matrix. In this case, the next two steps 
should be skipped). 
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c) Apply Chan's Pivoting scheme (steps 1 through 9 of the 
RROR algorithm) to II, Q, and R,, finding IL, Q, and R,. 


d) Let Il-IL, Q-Q; and R=R,. 


e) Identify the signal or noise subspace (Option to use a 
refinement as explained in Section B below). 


f) Apply the Minimum Norm to find the estimated source 
angles. 


g) Filter and store the source angles. 

The source code implemented ES generate the 
autocorrelation matrix is shown in Appendices K and L. The 
source code implemented to generate the adaptive algorithm is 
shown in Appendix D and the code corresponding to the rank-one 
modification shown in Appendix E. Appendix F presents the 
Minimum Norm algorithm. Finally, Appendix H presents the 
Minimum Norm identification procedure needed to isolate the 
Signal source locations. 

inesadaptive algorithm presented above is an alternative 
to the identification of the signal/noise subspaces via SVD or 
EVD decomposition. Note that steps 4 to 6 in Chan's algorithm 
are only needed when the element of maximum magnitude of the 
smallest right singular vector is not in the last position. 
Thus, additional reduction in computation time is obtained 
when no pivoting is needed. 

Next, we investigate how often pivoting is needed in 
Chan's algorithm. To test that effect, we ran two test cases 
with two sources (m=2). One of the sources is fixed, the 


other is time varying. Ten sensors are used to compute the 
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correlation matrix with 100 snapshots used to form the noise- 
free autocorrelation matrix and 200 updates are computed. The 
value of the SNR varies to try to identify correlation between 
the SNR and the need for pivoting during the RRQR decomposi- 
tion (step 4 of Chan's algorithm). In our test cases, if the 
Chan's pivoting scheme were applied blindly, steps 4 through 
9 would be executed for a total of 1600 times. (Because n- 
r+1=10-8+1=3 and the algorithm is iterated fromn to n-r+l, a 
total of eight times for each one of the updates is needed. 
Since we have 200 updates in our test case, we get 
8*x200=1600). 

First, we update Q and R directly using two successive 
rank-one modifications. Next, we update the nose james 
autocorrelation matrix using Equation (48) and perform a new 
complete RRQR decomposition (steps 0 through 9). Figures 5 
and 6 show the percentage of times that a pivoting is needed 
out of the 1600 tests as a function of the SNR (dB) of the 
Source. 

Figure 5 shows that the percentage of pivoting steps 
needed lies between 18% and 32% when two successive rank-one 
modifications are used to update the noise-free correlation 
matrix. This means that no pivoting is needed for every 
iteration. 

For the second test case, where a complete RRQR algorithm 
is applied in the updated noise-free autocorrelation matrix, 


the percentage of pivoting steps needed remains between 70% 
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PERCENTAGE OF PIVOTING NEEDED FOR 200 UPDATES 


PERCENTAGE OF PIVOTING NEEDED 
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Figure 5: Percentage of times that a pivoting is needed out of 
1600 iterations on the adaptive algorithm, using two rank-one 
modifications. 


PERCENTAGE OF PIVOTING NEEDED FOR 200 UPDATES 
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Figure 6: Percentage of times that pivoting is needed out of 
1600 iterations on the adaptive alg., updating the noise-free 
correlation matrix and applying a complete RRQR algorithm. 


and 889. The two successive rank-one modifications are more 
computationally expensive than the complete RRQR algorithm. 
On the other hand, they require less pivoting. In Section D 
below, we analyze the implications of this on the processing 


time for the two approaches. Furthermore, Figures 5 and 6 
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indicate that there is no correlation between the magnitude of 


the SNR and the need for pivoting. 


B. REFINEMENT ON THE RESULTS OF THE ADAPTIVE CASE 

This section presents an improvement to the signal 
subspace estimation procedure which is applied to the adaptive 
algorithm. The basic idea behind the improvement is based on 
the fact that the noise-free autocorrelation matrix, R,=QRII 
is Hermitian, and therefore R,-R,. So, R,{=IIR"'9". Recall that 
the signal and noise subspaces, Q, and Q, are contained in the 


matrix O. wPhuss 











Qj ON ir O 1 Of (51) 
H S S 
Rs OUR É) Os = DR qo Ms eso o)” Tax 
On 0 De R22 12 0 
and therefore 
H 
0 
11 
KO. y I, IL H 0 m II, R5 HL RS. (52) 


12 

The matrix resulting by the product R,Q, may be viewed as 
an one-step subspace iteration applied to the signal subspace 
Q,. This iteration scheme improves the results obtained as 
shown in the next section. The drawback is that the resulting 
iterated signal subspace R,Q, is no longer orthonormal. An 
additional orthonormalization step needs to be applied to the 
iterated signal subspace in order to use the Minimum Norm 
algorithm. Note that no reorthonormalization is needed when 


the MUSIC estimator is used. 
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C. SIMULATION RESULTS 

We now present the results generated by the simulation 
carried out for SNRs of 6, 10 and 20 dB. The noise is assumed 
to be zero-mean Gaussian and uncorrelated from sensor to 
sensor. We consider the case of two sources impinging on the 
ten-element array. The first source is assumed to be fixed at 
a normalized angle 0,-40?, the second source location @, is 
linear time-varying. Movement starts at 30° and stops at 
22.5°, after 100 snapshots are used to form the correlation 
matrix and 200 updates are used to Simulate the movement. 

Figures 7, 8 and 9 present the estimated DOA information 
obtained uSing the Eigenvector decomposition (EVD), the 
initial RRQR approximation, and the "refined" RRQR algorithm. 
The EVD decomposition is used for convenience instead of the 
SVD decomposition, as both span the same subspace. The 
results obtained for the refined RRQR technique are nearly 
identical to those obtained using the EVD technique. Table 25 
below shows means and standard deviations for the magnitude of 
the difference between the RRQR approximation for the signal 
DOA angle in degrees with/without refinement and the DOA 


generated by the EVD decomposition. 
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Figure 7: Position of the time-varying source for each update, 


SNR equal to 6 dB. 
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Table 25: MEAN AND STD. DEV. FOR THE MAGNITUDE OF THE DIFFER- 
ENCE BETWEEN THE RRQR APPROX. FOR THE SIGNAL DOA ANGLE IN DEG. 
FOR THE MOVING SOURCE WITH/WITHOUT REFINEMENT AND THE EVD. 


e W/O 
refinement 
(degrees) 


Average w 
refinement 
(degrees) 
tandard Devi- 
ation w/ re- 
finement (de- 
grees) 





The results are excellent even for the approximation without 
Boimement. The refinement improvement becomes better as the 
SNR increases. 

Results found for the largest principal angle between the 
projection of the signal subspace found via RROR and the true 
Signal subspace found via EVD are shown next. Figures 10, 11 
and 12 depict the results found for SNR values equal to 6, 10 


and 20 dB. 
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Figure 10: Largest principal angle between the signal subspace 


found via RRQR and the true signal subspace found via EVD for 
SNRz6 dB. 
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Figure 11: Largest principal angle between the signal subspace 
found via RRQR and the true signal subspace found via EVD for 
SNR=10 dB. 
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PRINC., ANGLE FOR THE SIG. SBSP-SNR=20 dB 


PRINCIPAL ANGLE IN DEGREES 
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Figure 12: Largest principal angle between the signal subspace 
found via RRQR and the true signal subspace found via EVD for 
SNR=20 dB. 


As can be seen, the results for the cases with refinement are 
much better than those without refinement. Table 26 shows the 
means and standard deviations obtained for the largest 
principal angle between the RRQR and EVD for both cases, with 


and without refinement. 
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Table 26: MEAN AND STANDARD DEVIATION FOR THE LARGEST 
PRINCIPAL ANGLE IN DEGREES BETWEEN THE RRQR APPROXIMATION AND 
THE EVD FOR THE SIGNAL SUBSPACE WITH/WITHOUT REFINEMENT. 


= SNR=6 dB [SNR=10 dB SNR=20 dB 


Average w/O a 3.4087 . 3442 
refinement 
(degrees) 

3 Ea 16502 "05965 
(degrees) E] 


Average W 2308 12252 T0023 
refinement 
(degrees) | sef 


tandard=Devi=s! 
ation w/ re- 
finement 
(degrees) 





Last, Figures 13, 14 and 15 present the behavior of the 
RROR approximation for the estimation of the DOA of the second 
Source that remained constant at 40°. Table 27 presents mean 
and standard deviation values for the magnitude of the 
difference between the DOA estimated by the RROR with/without 
refinement and the EVD. As can be seen, the RRQR results 


agree closely with EVD results. 
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Figure 13: DOA in degrees for the fixed source for each 
update, SNR equal to 6 dB. 
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Table 27: MEAN AND STD. DEV. FOR THE MAGNITUDE OF THE DIFFER- 
ENCE BETWEEN THE RRQR APPROX. FOR THE SIG. DOA ANGLE IN DEG. 
WITH/WITHOUT REFINEMENT AND THE EVD FOR THE FIXED SOURCE. 


9 YSNr=6 db [SNR=10 dB] SNR=20 GB — 
Average w/O re- 8970 3623 . 0418 
Standard Deviation .4603 .1619 90184 

w/o refinement 

Average w/ re- "0226 .0049 JO SUS ES 
Standard E T020 .0047 6.4206E-5 


/ refinement 
(degrees) 





















D. COMPUTATIONAL EFFICIENCY OF THE RRQR APPROXIMATION 

This section presents a basic estimation of the number of 
floating-point-operations (flops) needed to compute the RRQR 
approximation. The n-dimensional noise-free autocorrelation 
matrix is Square and is not considered complex valued for flop 
computation. This fact will be taken into account later. The 
number of flops necessary to perform one SVD decomposition is 
6n [REE ak Every time a new sample arrives, we need an 
O(n) operation to recompute the SVD [Ref. 13]. 

The RRQR algorithm is composed of three distinct parts. 
The computation of the initial QR factorization without 
pivoting, computed only once in the beginning of the algo- 
rithm; the computation of the least dominant right singular 
vector v of R, by inverse iteration, at each iteration; and 


the new QR factorization of R,,P, also at each iteration. 
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The first part takes 2n°/3 flops using the Householder 
algorithm if we do not need to accumulate Q. If Q is needed, 
as in our case, it takes 4n°/3 flops [Ref. 7]. The Modified 
Gram-Schmidt algorithm is preferred, since it takes only nº 
flops when the accumulation of Q is performed [Ref. 7]. 

Ignoring lower order terms, the second part of the RRQR 
factorization takes In’r flops to be iterated [Ref. 4], where 
r is the noise-free autocorrelation matrix rank-deficiency and 
I is the number of iterations used in the inverse iteration 
method. Our case uses I=3, therefore, the second part takes 
3n’r flops to be performed. 

In the third part of the RROR algorithm, 2n’r flops are 
needed when Givens rotations are used [Ref. 4]. However, note 
that not all elements below the main diagonal of the matrix 
R4,P needs to be annihilated because they are already zero. 
Therefore, a conditional "IF" statement may be used to verify 
if the element is already zero to save additional flops. 
Thus, the RRQR algorithm totals n°+5n’r flops at most. 

Following the first RRQR decomposition, we have two 
options. The first option updates the autocorrelation matrix, 
as in Equation (48) and performs a new QR decomposition. This 
update is followed by the pivoting scheme, as in steps 0 
through 9 of Chan’s algorithm. The second option updates the 
already obtained QR decomposition directly, using two 


successive rank-one modifications. A new Q, and R, is found 


p 


and the Chan’s pivoting scheme is applied. An analysis of the 
number of flops necessary to perform the two options is made. 

Recall that in the adaptive case, a double update must be 
performed for each sample. One update adds the new sample and 
the other subtracts the old one. If the rank-one modification 
option is used, the updates take 13n’ flops each for each 
snapshot, totaling 26n* flops [Ref. 7]. This does not include 
the number of flops necessary to proceed the pivoting scheme. 
If the complete RRQR factorization option is used, the QR 
decomposition using the Modified Gram-Schmidt algorithm takes 
n> flops. Note that in such a case, the autocorrelation 
matrix must be updated as R,,-R,4,*(x:xF),,,- (K° x") ag, which takes 


* for each multiplication and n^ for each 


4n" flops (n 
addition). The total is n'+4nº for each snapshot. 

Comparing n'+4nº with 26n’, we see that to update the 
noise-free autocorrelation matrix (finding R,, as in Equation 
(48)) and to take its QR decomposition is more economical than 
to perform two rank-one modifications for O « n « 22. 
Therefore, a complete RROR algorithm including a Modified 
Gram-Schmidt QR decomposition and the pivoting scheme is 
preferable when compared to the two rank-one modification 
updates, for a number of sensors smaller than 22. This method 
does not take into consideration the difference of pivoting 


schemes needed after the update. Table 28 presents a 


comparison for the number of flops needed for both processes. 
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Table 28: COMPARISON OF THE NUMBER OF FLOPS NEEDED TO THE 
SNAPSHOT UPDATE OF THE ADAPTIVE ALGORITHM VIA A COMPLETE RRQR 
FACTORIZATION AND TWO RANK-ONE MODIFICATIONS. 


[Snapshot update [Number of Snapshot update [Number of 
la complete flops need- [via two rank- flops need- 
RROR algorithm led one modif. ed 
| 4n First rank-one 
EMEN - Ain 
n Second rank-one 


Smallesterignt ade 
SV 
VRE ror É 20T | 


RiP 


Uerum Grizen 





As seen earlier, using two successive  rank-one 
modifications leads to pivoting for at most 32% of the 1600 
iterations used to perform the 200 updates. This is compared 
to the maximum E 88$ pivoting when updating the noise-free 
autocorrelation matrix and applying a complete RRQR algorithm. 
These numbers are used to compare the two approaches. 

For the third part of the RRQR algorithm, we need 2n’r 
flops. Using the first approach, two rank-one modifications 
over Q and R requires a total of 26n’+ (3+0.32x2)n’r flops for 
each update. In the second approach, the updating of the 
noise-free autocorrelation matrix and the application of a 
complete RROR algorithm requires n’+4n’+(3+0.88x2)n’r flops for 
each update. Figure 16 depicts the results showing the 
regions when the complete RRQR or the update approaches might 


be preferred. 
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ma at Two successive rank=one modificaktons mm 


¡is less expensive ; 


AUTOCORRELATION MATRIX SIZE (n) 


Campleté RRQR foctorizotidn 


is less expensive 





Q 5 10 VS 20 


RANK DEFICIENCY (r) 


Figure 16: Graph showing the regions where the RRQR or the two 
rank-one modifications approach is preferred. 


The MINNORM algorithm takes roughly a number of flops 
equals to 2nr plus the flops necessary to compute the 
polynomial roots. The root finding routine is iterative. 
Thus, it is impossible to obtain a specific expression to 
evaluate the number of flops necessary to run it. However, 
for the scenario simulated in the previous section, we were 
able to compute the mean and standard deviation of the number 
of flops spent by the MATLAB™ software to find the ten roots 
for each one of the 200 updates. 

To evaluate the sensitivity of the root finding routine to 
the SNR, we tested for SNR=-100, 6 and 100 dB. The results 


are shown in Table 29. 
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Table 29: MEAN AND STD. DEV. FOR THE NUMBER OF FLOPS NECESSARY 
TO FIND TEN ROOTS OF THE POLYNOMIAL FORMED BY THE MINIMUM NORM 
ALGORITHM. SNR=-100, 6, 100 DB. 





6 61007 2 BES 


64420 


According to Table 29, the number of flops necessary to 





the root finding routine presents a smaller standard deviation 
for SNR equals to 100 dB. This happens because at such a high 
SNR, the zeros are at more or less fixed locations. Thus, one 
can expect about the same number of iterations needed to 
identify them. The study of the polynomial root algorithm 
sensitiveness to the number of sources at different locations 
deserves further research and is out of the scope of this 
work. 

The REFINEMENT algorithm spends a total of 3n%+(1-2r)n?-rn 
flops, including the orthonormalization necessary to be used 
in conjunction with the MINNORM. When the problem is located 
under the line depicted in Figure 16, the updating via a 
complete RRQR algorithm is preferred in the most probable 
case.  Totalizing, the RRQR, the MINNORM and the REFINEMENT 


algorithms take a total number of flops of 4n*+(2.76r+5)n’+rn 


no 


for each update. This value was developed for a real valued 
noise-free autocorrelation matrix. To cope with future growth 
of the program, we might multiply the number of flops by a 4 
to 1 factor when implementing it operationally and by a 
roughly 4 to 1 factor to estimate the case of a complex 


matrix. 


E. INTEGRATING THE RRQR ALGORITHM INTO A REAL-TIME CASE 

Suppose we have a signal being received by a linear phased 
array on the earth surface from a moving object located at 10 
NM from the sensors at an angle of 80? from the sensors 


vertical (see Figure 17). 


m sources 


80° — 


e 
a 


Passive Linear 


e e 9 0€ .....e. o eoe e — Equispaced 
= h sensors Tay: 
m<n 


Figure 17: Situation of a signal being received by linear 
phased array sensors. 
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Also suppose that this signal is moving at an horizontal 
velocity v. Based on the simulations run in this thesis, the 
moving signal changed 5? in its DOA with respect to the array 
vertical. This was done in 200 snapshots. Each snapshot is 
taken at an interval of t, seconds. Therefore, the angular 


Speed of this moving signal is defined as 





E 
_ 1900m (53) 
POO). (ee ET 
The signal velocity is 
w'R 
ET 54) 
cos (0) 


where R is the distance between the signal and the array. 
Using (53) and (54) leads to: 


STR 


e EE 55 
= 36000'v'cos (0) un 


Assuming a signal is moving at sound speed as in Figure 18, 
0=80º, R=10 NM, a t, of 136.87 ms would be needed. 

For the scenario simulated in the previous section where 
n=10 and r=8, the number of flops needed for each update would 
be 16x6788=108.6 Kflops per update. This value is equivalent 
to 108,600/136.87 ms = 790 Kflops per second of microprocessor 
computing power. Assuming a reasonable computing power of one 
flop per clock pulse at 32 bits, a microprocessor would have 


to operate at approximately 790 KHz. Current commercial 
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microprocessors sold with a clock of 50 MHz would be able to 
implement this algorithm in real-time. Note that there is 
enough room to accommodate the root finding procedure 


neglected in the MINNORM algorithm flops computation. 


MS 


VII. CONCLUSIONS 

We have presented a fast algorithm to isolate signal and 
noise subspaces without performing the eigenvector 
decomposition of the autocorrelation matrix. This method is 
called the Rank-Revealing QR factorization. 

We have investigated the possibility of using the least 
dominant eigenvector instead of using the least dominant 
singular vector in step two of the RROR algorithm. 
Simulations have shown that the minimum singular vector 
Aa rated by the inverse iteration method gives better results 
than those obtained with the approximate smallest eigenvector 
generated by the same method. The Incremental Condition 
Estimator algorithm COn diag tan ao proxima loro the 
smallest singular vector has also been tested. Simulations 
have demonstrated that the inverse iteration again yields 
better results. 

The possibility of using a faster RRQR algorithm than the 
original algorithm with fewer iterations has been 
investigated. Simulations have shown that reducing the number 
of iterations worsens the signal/noise subspace estimations. 
Therefore, the original algorithm is preferred. 

Two spectral estimators have been tested to be used with 
the RROR algorithm, the MUSIC and the Minimum Norm. A limited 


number of Simulations have indicated that the Minimum Norm 
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resolved two signals for a lower SNR than the MUSIC method, in 
agreement with Rao [Ref. 12]. 

A relatively inexpensive computational refinement 
algorithm has been presented for the estimation of the RRQR 
Signal subspace. Simulations have shown that improvements at 
least as high as a factor of 20 are possible to be obtained 
for the signal angle of arrival, as compared with the original 
RRQR-based DOA results. Note that this refinement improvement 
becomes better as the SNR increases. 

An adaptive RRQR-based algorithm has been introduced to 
track the DOA of moving signals. Two options have been 
evaluated to compute correlation updates. The more adequate 
option may be determined depending upon the particular problem 
set up, e.g., the noise-free autocorrelation matrix rank 
deficiency (r) and the number of sensors (n). 

An evaluation of the number of flops required by the 
adaptive algorithm to find the DOA for two sources present has 
been determined. The results have shown the feasibility of 


the algorithm to solve a real-time problem. 
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Appendix A 


This function implements Chan's RRQR algorithm. 
Milton P. Ferreira, Sep/1992. 


moult parameters: 


r= matrix to which we want to apply the RROR factorization 
(noise-free autocorrelation matrix. 
a= rank deficiency of the matrix "r". 


Output parameters: 


Q= orthonormal matrix resulting from the RRQR factorization 
of Dae 

R= upper triangular matrix resulting from the RRQR 
REtorization of "r". 

e- permutation matrix (PI) resulting from the RRQR 
REcrorization of "r". 

nsbsp= noise subspace. 

ssbsp= signal subspace. 


function [Q,R,e,nsbsp,ssbsp]-rrgr(r,r1) 
n=size(r);n=n(1,1); 

ee-n-rl; 

wl-zeros (n); 

OER] =qr (r); 

e=eye(n); 

soun-ó;: 

for i=n:-1:n-r1+1; 


R11=R(1:1,1:1); % FIND THE THE LEADING ixi BLOCK (STEP 1) 
[u, sigmin,v]=ssvd(R11,3); % FIND THE LEAST DOMINANT 

FSINGULAR VECTORS AND SINGULAR VALUE (STEP 2) 
pt=eye (1); 


punfí-normív,inf); % FIND THE MAXIMUM ABSOL. VALUE 


$ ELEMENT OF THE LEAST DOMINANT SINGULAR VECTOR 
vauxl=abs (v) ; 
ind=find(vauxl==vinf):; 
End =1, % FIND THE POSITION OF THE MAX ELEMENT 
pues Otel Gia: % MAKE THE PERMUTATION TO PUT THE 
$ MAXIMUM ELEMENT AT THE i th POSITION (STEP 3) 
pe nd: )-ptl1, 2), 
pt (1, :) =vaux; 
wl(:,1)=[v;zeros(n-1,1)]; % ASSIGN v TO THE ith COLUMN 
SOFw (STEP 4) 
p-pt'; 
pud -zeros(i,(n-z)); 
peil? 1=zéros((n-1),1): 


s 


ptil22=eye((n-1) > 

ptil=[p ptili2 |= pe PN 

w-ptil'*w1; 

[Q1 Plikt Ser Re 

e=e*ptil; 

R12=R(1:1, 

if (n-1)== 
R22=[] 
else 
R22=R((1+1) :n, (1+1) :n); 

end; 

R= [R11til Q1'*R12 

OlI2szeros(1 n9) 

021=zeros((n-i),i) 

022-67 6 (Griese: 

0=0* [01 012 02 Cees. 


SEE x 
STEP 6 
STEP 7 


o? o? o? 


(Te) 0n); 
0 


. 
i 


zeros((n-i),i) R22]; 


. 
r 
. 
, 
. 
r 


end; 

end; 

nsbsp=Q(: ,n-rl+l:n): 
ssbsp=0 (27 tine 
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Appendix B 


wnet iom jik - givensl(x,y,n,i,k) 
%GIVENS Givens rotation matrix. 
% G = GIVENS (x,y) returns the complex Givens rotation 


C S X 


| 
| 
G = | 
| 
| 


| | 
| | 
| such that G * | 
| | 
| | 


cons) Vc Y 


il 
pa 


where c is real, s is complex, and c^2 + |s|^2 


Copyright (c) 1987-88 by The MathWorks, Inc. 
Modified by Milton P. Ferreira Sep/1992 


mput parameters: 


x= pivoting element [a(i,i)]. 

y= element we want to zero out [a(k,1)]. 

n= matrix dimension. 

i= column of the element we want to zero out. 
k= row of the element we want to zero out. 


Output parameters: 


jik= matrix "J". When multiplied by the matrix "a" will zero 
out the element a(i,k). 

absx = abs(x); 

fee absx == 0.0 

"5 00-953 To; 


o9 oV oV o o o oV o oV oV o o oV o ov 


a oe vi; 
c = absx/nrm; 

S - x/absx* (conj (y) /nzm) ; 
end 

jik=eye (n); 

soe (i, 1) =c; 

peek (kK, 1) =-conj (s); 

pei (1,k) =s; 

Elk (kK, k) =c; 
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Appendix C 


Computes Givens Rotations to zero out the elements of a 
matrix below its main diagonal. It is used at step 6 of 
the RRQR algorithm instead of a new QR factorization of 
R11*P. 

Milton P. Ferreira Sep/1992. 


Input parameters: 


w= matrix whose elements below the main diagonal we want to 
zero out. 


Output parameters: 


ql= new matrix Q, after applying Givens Rotations. 
rl= new matrix R, after applying Givens Rotations. 


FUN Ciel On ale guga) 


1=size(w) ;1=1(1,1); 
qt=eye (i); 
jo pale (12 31 
for p=l mingled spe. 
if w(q,p)”=0, 
jpq=givensl (wp på vig PTE 
qt-jpa*qt; 
w-Jpq*w; 
end; 
end; 
end; 
qu qu 
ri=-w; 
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Appendix D 


Algorithm for the adaptive tracking of a moving source 
Main program. 
Milton P. Ferreira Sep/1992. 


miput parameters: 
db= SNR. 
OWEPUt parameters: 


RO= autocorrelation matrix. 

Rl= noise-free autocorrelation matrix. 

Dach column of Y is one output vector "x" from the "n" 
Sensors. 

nest- # of snapshots used to compute the autocorrelation 
matrix. 

nupd= # of updates used to simulate the movement. 

ipp- autocorrelation matrix dimension. 

nsin= number of sources. 


o? opp opp oo opp opp opp opp opp op o? o? o? opp opp oe opp a? o? o? 


ROB RI,Y,nest,nupd, ipp,nsin]=cor4ámi (db); 
[Q,R,e,nsbsp,ssbsp]=rrar(Rl,ipp-nsin); 
for i=nest+l:nest+nupd; 
old=Y(: ,i-nest): 
OR] =givgr OE, CMM ord) 
new=Y(:,1); 
lO, RX] =givgr (0O,R,e,new, new) ; 
[Q,R,e,nsbsp,ssbsp]=rrqe(Q,R,e,nsin,ipp); % is the RROR 
% subroutine without the initial QR decomposition 
[mags,angs] -minn(nsbsp,ssbsp,ipp); 
[ma,an]-ident (mags,angs); 
angupd (i-nest)-an(2,1); 


end; 

mrot(1,angupd) ; 

End; 

title(['SIGNAL SUBSPACE - SNR-',int2str(db),' dB']); 


ylabel('SOURCE LOCATION'); 
xlabel('NB OF UPDATES'); 
le 7, .5,' SOLID =UPDATE','sc'); 
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Appendix E 


Computes the rank-one modification of a square matrix. 
Tnput parameters: 


Q= orthonormal matrix resulting from the RRQR 
factorization? 

R= upper triangular matrix resulting from the RROR 
factorization: 

e= permutation matrix (PI), resulting fron the RROR 
factorization. 

u and v= u and v vectors that modify the matrix to be 
updated. 


Output parameters: 


Q1= matrix Q after modifying. 
Rl= matrix R after modifying. 


Milton P. Ferreira Sep. 1992. 


o? oV oV o o o? AO AL A o? o A A A o? o? oV A o a 


function [DINPRIS GI VGA OR ENE 

w=Q*u; 

Q1-0; 

H=R; 

i-size(w);isi(1,1); 

for ge uq 
Jpa=givensiiw(1,1),w(q, 
w-jpq*w; 
H-jpq*H; 
Q1-Q1*jpq'; 

end; 

R1=H+w*v'*e; 

for qd su 
jJpa=givensi (Ri(la, a) RITO M IUE 
Ri=jpg*R1; 
Q1=Q1*jpq'; 

end; 
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Appendix F 


function [mags, angs] -minn (nsbsp, ssbsp, ipp) 


oV oV o o? o o? o o o oV o oV o o oV o 


compute the noise and signal zero locations using the TK 
min-norm 
alg. (version 1.0 10/14/91), Monique P. Fargues. 


Muu: parameters: 

ipp= correlation matrix dimension 
nsbsp= noise subspace. 

ssbsp= signal subspace. 


Output parameters: 


mags= magnitude of the polynomial roots. 
angs= phase of the polynomial roots. 


@m=Zero0s(1:ipp) ;ds=zeros(1l:ipp); 
clear g 

SuEnsbsp(í1,:); $noise zeros 
En=nsbsp (2: KR Å 

sie ipp) =En*g'/(g*g'); 

anml) =l; 

clear g 

g=ssbsp(1,:); %signal zeros 
Es=ssbsp(2:ipp,:); 

ESI/(1-g*g'); 

@s(2:ipp) =Kg* (Es*g") ; 

Gent) =1; 

flops (0) ; 

anr=roots(dn); 

dsr=roots (ds); 

mags=abs (dsr) ; 

angs=angle (dsr) *180/pi; 
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Appendix G 


Computes the "refined" signal subspace 
Version 1.0 08/21/92, Monique P. Fargues. 


Input parameters 


R= upper triangular matrix resulting from the RRQR 
tactorizataome 

e= permutation matrix (PI) resulting from the RRQR 
factorization. 

nsin= # of sources. 

ipp= correlation matrix dimension. 


Output parameter. 


ref= "refined" and reorthonormalized signal subspace. 


function ref=refine(R,e,nsin, ipp) 
pil=e(: 1- ASIn) 

pi2=e(: nsinti: Ipp); 
RIISR(lI:nsgymnou Ms E 

R12=R (1:nsin,nsin+1:ipp); 
ref=orth(pil*R11'+pi2*R12'); 
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Appendix H 


identify the signal magnitude and location of the roots 
version 1.0 10/13/91, Monique P. Fargues. 
UE parameters: 


mag= magnitude of the polynomial roots 
ang= phase of the polynomial roots 


Output parameters: 
mag_Y= vector containing the magnitude of the roots 


corresponding to the sources. 
ang r= idem to the phases. 


A AL AL AL A o? o? o? o? o? o? o o? 


function [mag r,ang r]=ident (mag,ang) 
[m,n]=size(ang); 

ses ds 

ERRO (1:2,1)=[0;0];mag r0(1:2,1)=[0;0]; 
mene 1=1:m 

if mag(i) g(i)<50 & ang(i)>10 %look at both 
DuC 
1): 


if k==3, break, end 

end 

$ sort to insure proper separation of sines 
Pmeger(:,1),lJ=sort (ang r0(:,1)); 
mager(:,1)=(mag rO(1,1)); 
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Appendix I 


function [sigmin, x] =icest (R) 

1/7/92 **kkkk-icest.m-******x 

version 1.0, Monique P. Fargues. 

incremental condition estimator (modified from Bischoff 
paper) 

computes an estimate for minimum singular value and 
Vector 

associated 

with an upper triangular matrix 

initial matrix | gamma v 

| 0 R | 


R= matrix we want to estimate the smallest S. value and 
its S. vectors. 


output parameters: 


sigmin: minimum singular value 
Xx: IS OS singular vector 


% 
$ 
$ 
$ 
% 
% 
$ 
$ 
$ 
$ 
* Input parameters: 
% 
% 
% 
% 
$ 
% 
% 
6 
o 


[m,n] =size(R) ; 

$if(m”=n), error ('R is not square'), end 

sif (any (diag(R)==0)) % matrix is singular 

% smin=0: vmin= zeros (n,D rvmin(un CE G eal 
% return 

send 


x(n,1)=1/R(n,n) ; 


for icn- ren 
xx(L:n-1 =x (irl mI 
v=R(i,i+1:n)'; gamma= 
alpha=v' *xx; 

if alpha”=0 
beta= (abs (gamma) *norm(xx,2))*2 + abs(alpha)*2 -1; 
eta=beta/ (2*abs (alpha) ) ; 
sqr=sqrt (eta*2+1); 
nu=eta + sqr ; % sart(eta”2+1); 
temp=abs (alpha) *nu; 
root=temp+l; sqr-sqrt (nu^2-*1); 


JE 
RJ) 


s=temp/(alpha*sar) : * sqrt(nu^241)); 
c--1/sqr ; % sqrt (nu*2+1); 
else 
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root= (abs (gamma) *norm(xx,2)) "2; 
if (root>1), c=1, s=0; 
else, c=0, s=1; end 

end 


x(i:n,1)=[(c-stalpha) /gamma;s*x(1+1:n,1)]; 
end mor initial loop 
Euorm-normí(xí(i:n,1),2); 


seen 1) =x (1:n,1) /xnorm: 
Sromin=1/xnorm; 


gt 


Appendix | 


function [vsv, sigmin,usv]=ssvd(a,k) 


THIS FUNCTION IS THE IMPLEMENTATION OF THE INVERSE 
ITERATION 

TECHNIQUE TO FIND THE SINGULAR VECTORS AS SHOWN IN THE 
PAPER 

"DEFLATED DECOMPOSITION OF SOLUTIONS OF NEARLY SINGULAR 
SYSTEMS" - TONY F. CHAN - PG 746 - SIAM J. NUMER. ANAL. 
VOL 21 #4 AUG. 1984 

[vsv, sigmin,usv]=ssvd(a,k). 


Input parameters: 

a= matrix we are looking for the smallest singular 
Net OE 

k- number of iterations desired. 


Output parameters: 


usv, vSv= right and left smallest S. vectors. 
sigmin= minimum singular value. 


Milton P. Ferreira 


oV? o? o9 o9? AS AL AL AL A A o? o? o o o? o o? o? o? o? o? o9 


n-size(a); 

n=n(1,1); 

v=ones (n,1); 

for is1:k; 
util=a\v; 
u=util/norm (util); 
vtil=a'\u; 
vevtil/norm(vtil) ; 

end; 

VSV=V; 

util=a\vsv; 

usv=util/norm(util): 

sigmin=1/norm(util); 
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Appendix K 


Mmmetion {RO,R1,Y,nest,nupd, ipp,nsin) =cor4mi (db) ; 
Monique P. Fargues. 


COMPUTE THE CORRELATION FUNCTION ONLY 
Er rocplitz correlation function 


EuDUt parameters: 
øp= SNR. 
puput parameters: 


RO= autocorrelation matrix. 

Rl= noise-free autocorrelation matrix. 

REI cn column of Y is one output vector "x" from the "n" 
sensors. 

nest- # of snapshots used to compute the autocorrelation 
matrix. 

nupd= # of updates used to simulate the movement. 

ipp= autocorrelation matrix dimension. 

nsin= number of sources. 


A AY AL AL A o oV oV oV oV oV o AO AD AO AO AO Al A a a o 


ENG format compact 
seed1=1042;rand('seed',seed1) 


GO r=0 ; Fampue (rue est correl 1/0: '); 
itop=0; camp pisos plz mon toeplitz 1/0209!) 
nupd=200; sinput ('update nb nupd: *); 


$nupdO-input('update nb nupdO (drop in angle): '); 
EXE -q-5; % input ('del freq (in percen*freq(1): '); 
COL == 

Mortimer true cor. seq\n’') 


Porina iest. cor. seqin') 


Mest=100'; 
eeb-osnput(' input gdb: {x x] '); 
gdb= (dD, db] 1 
ang= [30 40]; 
freqt-[7.5 9]; 


Source angles 
temporal frequencies 


o? oV oV o? o? 


Nor. number of sensors 
nsin=2; number of sources 
ssig=1.; ref. noise varlance 
ne sarte (-1): 

sigma=1; 
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9 


for MS: nevn % linear amplitude 

A(i)=sart (2.)* (10 (PN ces 
end 
EEE 
% compute the true correlation sequence 
for kl pE 
res=0; 
if k== 

res-ssig^2; 
end 

FOr =i Noam 

res=res+ (A(i)%*2) *exp(jc*(k-1) *ang(i) *pi/180.)/2.; 

end 
R (k) 2res; 
end 
RO-toeplitz(R); 
else 


compute the estimated correlation seq 
based on nest data points 
rand ('normal') 
%sigm=sqrt(10"((sigma)/10)) 
For i= osin 
A(1i)=sqrt (2) *sigma* (10” (gdb(i)/20)) 
end 
Eon pE 


o? oo 


Y(1i,1:nest+nupd) =sigma* (rand(1,nest+nupd) +jc*rand(1,nest+nup 
ae: 


end 

freq-ang*pi/180; 

rand( (uniform) % create uniform variable dist. 
sn (ppa) 

X1=pi* (rand(1,nsin*ipp* (nest+nupd) )-.5); 

FreqQ=freq (1); 

for j2=1:nest+nupd % j2: time 


snapshot 

sfreq(1)=freqO - (pi/180)*(j2-1)*del freq/nupd; $del freq 
% change in nupd samples 

+1f j2>nest+nupd0, freq(l)=freq0*del freq; end *step freq. 


% change 


number of sines 


for i-1:nsin e 
ji: sensor position 
* (1 


FOr JIS DO 


. o9 o? 


temp=(j1*freq(1)+j2*freqt (1) +X1(1 sj2* (1-1) app jae 
Y(91L,72)-Y (IL, Nl | *exp (jcttemp) ; 

end 

end 


end 
RO=Y (271: mest) *Y 0 a mese 
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end 


get the noise free version of RT & normalize the matrix 
compute the EVD of the noisy matrix for later computations 
Rom original correlation function 
[Ve, De] =eig(RO) ; 
sort the eigenvalues 
s (Des, 1] =sort (diag (De) ) ;12=flipud(l) ; 
$Ves=Ve(:,12(:)); sorting in descending order 
$Vnoi=Ves(:,nsin+1:ipp); isolate the noise vectors 
+Vsig=Ves(:,1:nsin); | | | &$ ----------- signal ------- 
R1-R0-nest*ssig^2*eye (R0); noise free correlation 
matrix 
$R1=R1./R1(1,1); % potential matrix 


o 


$ normalization 


V oV o A o? o? 


A A o o? 


BNUcor--1, 


R1-R0-ssig^2*eye(R0); % noise free correlation 
matrix 
end; 
mei cor==0, 

R1-R0-nest*ssig^2*eye (RO); 
end; 
ment op==1, 

Rl=toep(R1); 

end 
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Appendix L 


Transforms a matrix to a Toeplitz matrix. Used with 
"cor4mi.m". 


Tnpuüut parameter: 

Matrix to be transformed. 

JQUEPUEA parameters 

sum= Toeplitz matrix after transformation. 


Milton P. Ferreira Sep/1992. 


oV oV o? o? o? o? o? o? o? o? o o? o? 


function sum-toep (R1); 

n=8128 (RI nn IOE 

sum-zeros (R1); 

for i=-n+1:n-1; 
a=diag(R1,1); 
s-mean (a); 
sl=size(a); 
SIsseq p DF; 
s2=ones (s1,1)*s; 
sum=sum+diag(s2,1); 

end; 
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Appendix M 


Computes the smallest eigenvector by inverse iteration. 


mole parameters: 

a= matrix from which we want to compute an approximation 
of the smallest eigenvector. 

u= # of iterations. 


Output parameter: 


X= approximation for the smallest eigenvector. 


o? o? o? o? o? o? oV o? oV AP A AP DL do 


Milton P. Ferreira Sep/1992. 


mEDCLIOn x=eeig(a,u) 
nsesize (a); 
n-n(1,1); 
ones (n, 1); 
Gor i=l:u; 
t=a\x; 
x=t/norm(t,inf); 
end; 
x=x/norm (x); 
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